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引言 

首先，我们介绍一下使用计算机进行数值计算的几个基本术语。众所周知，计算 
机使用2进制表示数字，通常分为32位 （bit) 和64位两种。如果是32位，有一位表 
示符号，7位用来记录方幂，剩下24位用于分数，这样一来，每个计算机表达的数和 
原先给的数之间的精度为2的负25次方，大约是8位十进制小数，这就是常说的“单 
精度”。如果是64位，有一位表示符号，11位用来记录方幂，剩下52位用于分数，这 
样以来，精度为2的负53次方，大约是16位十进制小数，这就是“双精度”。所谓浮 
点运算，都是将每一个参与计算的数字表示成32或64位的浮点数来进行的。不难看出， 
使用单精度和双精度，计算过程产生的舍入误差是不一样的。 

其次，读者还需了解目前广泛使用的两个商业数学和工程计算 软件： MATLAB 和 
Maple。MATLAB 是矩阵实验室 （Matrix Laboratory) 的简称，由美国 Math Works 公 
司出品，它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿 
真等诸多强大功能集成在一个易于使用的视窗环境中，在工业设计、科学计算以及课 
堂教学等方面都有非常广泛的应用。相比于 MATLAB，Maple 的特点是符号运算，其 
功能覆盖诸多的数学领域，如微积分、线性代数、微分方程、积分变换、概率论和数 
理统计、图论、张量分析、微分几何、线性规划、组合数学、抽象代数、泛函分析、 
数论、复分析和实分析、特殊函数、编码和密码理论、优化等。与历史更久的 Fortran 
和 C 等计算机语言相比， MATLAB 和 Maple 是更高级、更加人性化的语言，其语法相 
对简单，经常是直接输入数学公式即可达到计算目的。 

有了这些准备，现在言归正传。对下面这个定积分 
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第一台可编程计算机 Z3, 包括浮点运算（在慕尼黑的德意志博物馆展出的复制品） 

分步积分，我们可以得到以下的递推公式 

I n = l~nl n _ l , "二 1 , 2 , ... 

这里 / Q 二1 — q ~ 1 o Renata Babuskova 1 早在1964年就曾讨论这个迭代的数值稳定性。 
以单精度计算，几步递推之后就出现了负值。在今天的标准 PC 上用 MATLAB 计算（双 
精度），我们观测到最初的一个递减序列直到/ 16 二0.0374。然而，随后的结果却令人 
惊愕： 


/ 17 二 0.3259, / 18 = -5.1930, 1 19 — 104.8628, …， 

我们看到一个迅速发散的交错序列。究其原因，最初的每一步，两个非常相近的数相减， 
只有很少几位有效 数字； 而这个数值误差又被其后的迭代以每步《倍的速度放大，从 
而产生了一个迅速发散的序列。 

如果我们改变策略，用以下快速收敛的级数计算 

7 = J _1_ + _1_ 

n ~ ^TT ~ {n + l)(n + 2) + 0 + 1)0 + 2 )(n + 3) — ••‘ 

则可以得到任意精度的逼近值。 

类似的例子不胜枚举，我们可以在1966年的一本书中 2 找到下面的 例子： 

... (((((1+2) .2) + 3) • 3) + 4) . 4… 

不断进行先除后乘的运算。 Karel Segeth 在不同的计算机上就这个例子进行了成千上万 
次乘除法运算，得到了许多有趣的结果。 

基于这样一个开场白，此文通过几个精选的例子警醒世人，在解读计算机得到的 
数值结果时，须要格外谨慎。 

一个递减序列会被浮点运算变为递增序列吗？ 


1 R. Babuskova: Uber numerische Stabilitat einiger Rekursionsfor mein. Apl. Mat. 9 (1964), 186—193. 

2 I. Babuska, M. Prager, E. Vitasek: Numerical Processes in Differential Equations. John Wiley & Sons, 
London, New York, 1966. 
























Maple 的标志 MATLAB 的标志 

取 a Q = 1, 屮 = 11 — S 开始以下递推： 

34 3 

a n+2 = ~[i an+l ~Yi an， ^ = 0,1,2, ••- 

准确解％二 111 是一个单调递减序列。然而，使用 MATLAB 在标准 PC 上计算， 
我们只能看到递减序列到 a 12 二 0.2068 *10— 11 为止，随之序列开始递增，到 a 4Q 二 
40.44,然后迅速增长。这个例子说明即便很小的两个相近的数相减也会带来大麻烦。 
有趣的是，用以下稍加变化的递推公式， 

a n+ 2 = ^-( 34 ^ + i - H)，w = 0,1 ， 2, ... 

再用 MATLAB 计算，可以看到序列从 a 12 开始变为负值，但绝对值迅速增长。 

这里我们给出序列中第50个数在不同计算精度下 的值： 
a 5( ^2*10 n 单精度，序列在《二8变为递增 序列； 

a 50 ^10 10 实值精度（相当于12位十进制精度），序列在《二9开始 递增； 
a 50 ^5*10 5 双精度，序列在《二12开始 递增； 

fl 5 0~ !° 3 扩展精度（相当于20位十进制精度），序列在 W 二14开始递增。 

Babuska 的例子 

考虑定义在 [0, 1] 区间上的狄利克雷函数/:在有理点/二0,而在无理点/二1。 
由于这个区间上所有无理数组成的集合的勒贝格测度为1，/的勒贝格积分 

1 

^ f(x)dx = l, 

0 

然而采用任何精度的数值积分都将得到 

i=l 

因为由浮点运算产生的 A 值全部为有理数。 

这个例子比较极端，数值积分的误差永远为1，原因在于/不是一个光滑函数。而 
数值积分收敛理论总是要求被积函数/具有一定的光滑性。 
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Muller 的例子 3 

取 q 二 e - 1 二 1.718281828 - 产生一个表面看来简单而又平常的 序列： 
a n = n ( a n _ x - 1), « 二2, 3,… （1) 

用数学归纳法可以 得到： 

/ 1 , 1 . 1 .) 
n Oz + 1)! Oz + 2)! J 

显然 { an } 是一个极限为1的递减序列。使用 MATLAB 数值实现（1)，开始我们的确看 
到递减，但好景不长，从 a 15 = 1.0668 起，奇怪的事发 生了： 

a l6 — 1.0685, a l7 — 1.1652, a ls = 2.9731, a l9 — 37.4882, a 20 — 729.7637, ... 

以后的序列数迅速膨胀。 

下面我们记录了序列中第25号元素在不同运算精度的 数值： 

a 25 ~- 1.306946321 • 10 13 实值精度， 

« 25 ~ + 1.201807248 • 10 9 双精度， 
a 25 ~- 7.3557319606 • 10 6 扩展精度。 

现在我们来看看提高计算有效位数^ (十进制）的 效果： 

D a 25( 乃） 

20 615990.413139 

21 -4457.98859386 

22 -4457.98859386 

23 195.374419140 

24 40.2623187072 

25 -6.27131142281 

26 1.48429359885 

注意，对应于 D 二21和 Z ) 二22的值相等，原因是欧拉数 e 的第21位是0。当 Z ) 增 
加到26以后，我们得到的计算值开始“靠谱” 一一接近真值 a 25 二 1.03993872967 … 
比如 a 25 (30) 二 1.039897 … 

让我们进一步探 if 一下这一奇怪现象的原因。以~代表第/步的截断误差，得到 
心二 e — 1 + q 二 q 

= 2 { a x _ 1) + = 2 { a x + q — 1) + e 2 = a 2 + 2!^! + e 2 , 


«25 z 25(« 2 4 - 1) + e 25 = a 25 + 25!q + e 3 + ... + e 25 • 

因而，累计截断误差 

a 25 - 及25 = _25 !〔^ + + ^ + ... + j ， 

它的大小特别依赖于最初的几个截断误差 A , 4 ，…这个例子再一次说明在数值计算中 
两个非常接近的数相减是多么的危险。 

3 J.-M. Muller et al.: Handbook of Floating-Point Arithmetic. Birkhauser, 2009. 












